#include "graph.h"
#include <QtCore>
#include <QtGui>
#include <QVector>
#include <qwt_scale_draw.h>
#include <qstringlist.h>
#include <qfile.h>
#include <qtextstream.h>
#include <qwt_scale_widget.h>
#include <limits>
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
#include <cmath>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <utility>
#include <algorithm>

#include "Numerical_Pde.h"
#include "Numerical_Optimization.h"
#include "Hydraulics.h"
#include "VestaData.h"
#include "oilpipe.h"
#include "gaspipe.h"
#include "gravity_simple.h"
#include "pig.h"

using std::vector;
using std::map;
using std::ofstream;
using std::pair;
using std::copy;
//using namespace Eigen::Sparse;

using namespace pdelib::force_flow;
using namespace pdelib::gravity_flow;
using namespace pdelib::combined_flow;
using namespace pdelib::gas;
using namespace pdelib::force_nonisotherm;
using namespace pdelib;

graph::graph(QWidget *parent, Qt::WFlags flags)
	: QMainWindow(parent, flags)
{
    ui.setupUi(this);
    plot.reset(new SimpleQwt(ui.Plot));

	connect(ui.btnSolver, SIGNAL(clicked), this, SLOT(on_btnSolver_clicked));
	connect(ui.slTime, SIGNAL(value_changed), this, SLOT(on_slTime_valueChanged));
}

void graph::get_setting_from_GUI(oilpipe_parameters* params, 
    pressure_series_equation* p_left, pressure_series_equation* p_right)
{
    //double p0 = ui.sbP0->value()*1.013e5 + M_G * params.density * params.profile.height_begin();
    //double pl = ui.sbP1->value()*1.013e5 + M_G * params.density * params.profile.height_end();

    /*double p0 = ui.sbP0->value() * M_G * params.density;
    double pl = ui.sbP1->value() * M_G * params.density;

    vector<double> dt(3), p(4);
    //double dp = ui.sbDP0->value()*1.013e5;
    double dp = ui.sbDP0->value() * M_G * params.density;
    p[0] = p0;  p[1] = p0 + dp; p[2] = p0 + dp; p[3] = p0;

    dt[0] = 20; dt[1] = 200;    dt[2] = 20;
    pressure_series_equation left_force_boundary(p, dt);
    //pressure_series_equation right_force_boundary(pl, pl + ui.sbDP1->value()*1.013e5, 20);
    pressure_series_equation right_force_boundary(pl, pl, 20);*/

}

void graph::initGraph(){
	//ui.Plot->setTitle(QString::fromLocal8Bit("Pressure"));
    //legend
}


graph::~graph()
{ 
	//delete ui;
    /*delete leg;
    delete grid;
    delete curv1;
    delete curv2;*/
}


struct ode_system {
    virtual MatrixXd system_matrix(double x, const VectorXd& U) const = 0;
    virtual VectorXd dissipation_function(double x, const VectorXd& U) const = 0;
    virtual unsigned equation_count() const = 0;
};

struct ode_gaspipe: public ode_system, public pipe_parameters, public gas_properties {
    ode_gaspipe(double _massflow, pipe_parameters pipe, gas_properties gas, EoS_zfactor& _eos) 
        : pipe_parameters(pipe)
        , gas_properties(gas)
        , massflow(_massflow)
        , eos(_eos)
    { }

    double coeff_a1 (double P, double T) const {
        double S = circle_area(diameter);
        double a1 = 1/sqr(massflow/S) - (M_R * T / sqr(P)) * (eos(P, T) - P * eos.derivative_pressure(P,T));
        return a1;
    }

    double coeff_a2 (double P, double T) const {
        double a2 = -M_R*T/P  * (eos(P,T) + T * eos.derivative_temperature(P,T));
        return a2;
    }

    double coeff_b1 (double P, double T) const {
        double b1 = M_R/P * (eos(P,T) + T * eos.derivative_temperature(P,T));
        return b1;
    }

    double coeff_b2 (double P, double T) const {
        double b2 = heat_capacity_pressure;
        return b2;
    }

    double coeff_c1 (double P, double T, double x) const {
        int index = x / profile.segment_length;

        double S = circle_area(diameter);
        double density_working = eos.get_density(P, T);
        double velocity = massflow / (density_working * S);

        double dynamic_viscosity = gas_dynamic_viscosity_NTP(get_reduced_pressure(P), get_reduced_temperature(T));
        double re = density_working * velocity * diameter / dynamic_viscosity;
        double lambda = hydraulic_resistance_mitichkin(re, relative_roughness);

        double c1 = -1/sqr(massflow/S) * (
            + 0.5*lambda*sqr(massflow/S) / (diameter*density_working) 
            + density_working * M_G * profile.sinus[index]
        );

        return c1;
    }

    double coeff_c2 (double P, double T) const {
        double S = circle_area(diameter);
        double density_working = eos.get_density(P, T);
        double velocity = massflow / (density_working * S);
        double dynamic_viscosity = gas_dynamic_viscosity_NTP(get_reduced_pressure(P), get_reduced_temperature(T));
        double re = density_working * velocity * diameter / dynamic_viscosity;
        double lambda = hydraulic_resistance_mitichkin(re, relative_roughness);
        return -M_PI * diameter * ambient_heat_transfer*(T - ambient_temperature) / massflow 
            +0.5*lambda*sqr(massflow/S) / (diameter*sqr(density_working));
    }

    virtual MatrixXd system_matrix(double x, const VectorXd& U) const {
        double P = U(0), T = U(1);
        MatrixXd A(2,2);
        A(0,0) = coeff_a1(P, T);
        A(0,1) = coeff_b1(P, T);
        A(1,0) = coeff_a2(P, T);
        A(1,1) = coeff_b2(P, T);
        return A;
    }
    virtual VectorXd dissipation_function(double x, const VectorXd& U) const {
        double P = U(0), T = U(1);
        VectorXd b(2);
        b(0) = coeff_c1(P, T, x);
        b(1) = coeff_c2(P, T);
        return b;
    }
    virtual unsigned equation_count() const {
        return 2;
    }
    double massflow;
    EoS_zfactor& eos;
};


struct ode_gaspipe_isotherm: public ode_system, public pipe_parameters, public gas_properties {
    ode_gaspipe_isotherm(double _massflow, pipe_parameters pipe, gas_properties gas, EoS_zfactor& _eos) 
        : pipe_parameters(pipe)
        , gas_properties(gas)
        , massflow(_massflow)
        , eos(_eos)
    { }

    double coeff_a1 (double P, double T) const {
        double S = circle_area(diameter);
        double a1 = 1/sqr(massflow/S) - (M_R * T / sqr(P)) * (eos(P, T) - P * eos.derivative_pressure(P,T));
        return a1;
    }

    double coeff_a2 (double P, double T) const {
        double a2 = -M_R*T/P  * (eos(P,T) + T * eos.derivative_temperature(P,T));
        return a2;
    }

    double coeff_b1 (double P, double T) const {
        double b1 = M_R/P * (eos(P,T) + T * eos.derivative_temperature(P,T));
        return b1;
    }

    double coeff_b2 (double P, double T) const {
        double b2 = heat_capacity_pressure;
        return b2;
    }

    double coeff_c1 (double P, double T, double x) const {
        int index = x / profile.segment_length;

        double S = circle_area(diameter);
        double density_working = eos.get_density(P, T);
        double velocity = massflow / (density_working * S);

        double dynamic_viscosity = gas_dynamic_viscosity_NTP(get_reduced_pressure(P), get_reduced_temperature(T));
        double re = density_working * velocity * diameter / dynamic_viscosity;
        double lambda = hydraulic_resistance_mitichkin(re, relative_roughness);

        double c1 = -1/sqr(massflow/S) * (
            + 0.5*lambda*sqr(massflow/S) / (diameter*density_working) 
            + density_working * M_G * profile.sinus[index]
        );

        return c1;
    }

    double coeff_c2 (double P, double T) const {
        double S = circle_area(diameter);
        double density_working = eos.get_density(P, T);
        double velocity = massflow / (density_working * S);
        double dynamic_viscosity = gas_dynamic_viscosity_NTP(get_reduced_pressure(P), get_reduced_temperature(T));
        double re = density_working * velocity * diameter / dynamic_viscosity;
        double lambda = hydraulic_resistance_mitichkin(re, relative_roughness);
        return -M_PI * diameter * ambient_heat_transfer*(T - ambient_temperature) / massflow 
            +0.5*lambda*sqr(massflow/S) / (diameter*sqr(density_working));
    }

    virtual MatrixXd system_matrix(double x, const VectorXd& U) const {
        double P = U(0), T = U(1);
        MatrixXd A(2,2);
        A(0,0) = coeff_a1(P, T);
        A(0,1) = coeff_b1(P, T);
        A(1,0) = coeff_a2(P, T);
        A(1,1) = coeff_b2(P, T);
        return A;
    }
    virtual VectorXd dissipation_function(double x, const VectorXd& U) const {
        double P = U(0), T = U(1);
        VectorXd b(2);
        b(0) = coeff_c1(P, T, x);
        b(1) = coeff_c2(P, T);
        return b;
    }
    virtual unsigned equation_count() const {
        return 2;
    }
    double massflow;
    EoS_zfactor& eos;
};


class ode_solver_euler {
public:
    ode_solver_euler(ode_system* _ode) 
        : ode(_ode)
    { }
    void solve(const VectorXd& initial_condition, VectorXd* buffer_begin, VectorXd* buffer_end, double step_size) {
        unsigned segment_count = buffer_end - buffer_begin;
        buffer_begin[0] = initial_condition;
        for (unsigned index = 1; index < segment_count; ++index) {
            MatrixXd A = ode->system_matrix(step_size * (index - 1), buffer_begin[index - 1]);
            VectorXd b = ode->dissipation_function(step_size * (index - 1), buffer_begin[index - 1]);
            buffer_begin[index] = buffer_begin[index - 1] - step_size * A.inverse() * b;
        }
    }
private:
    ode_system* ode;
};


void test_gas_stationary() {
   
}


void graph::on_btnSolver_clicked() {
    test_gas_stationary();

    return;

    pipe_parameters pipe;
    pipe.diameter = ui.sbDiameter->value();
    pipe.ambient_heat_transfer = 0;
    pipe.ambient_temperature = 0;
    pipe.relative_roughness = 0.5e-3;
    pipe.wall_thickness = 5e-3;
    
    pipe_profile = profile_2_points(ui.sbWellLength->value(), ui.sbWellLength->value());
    pipe.profile =  subdivide_profile(heights_to_angles(pipe_profile), ui.sbSectionCount->value());

    gas_properties gas; // default natural gas
    EoS_ONTP eos(gas);
    ode_gaspipe model(ui.sbMassFlow->value(), pipe, gas, eos);

    ode_solver_euler solver(&model);
    
    VectorXd initial_U(2);
    initial_U << ui.sbP0->value() * 1e5, ui.sbT0->value();
    QVector<VectorXd> buffer(ui.sbSectionCount->value(), VectorXd::Zero(2));

    solver.solve(initial_U, &buffer.front(), &buffer.back() + 1, pipe.profile.segment_length);
	
    //ui.Plot->setAxisScale(QwtPlot::xBottom, 0, buffer.size() * pipe.profile.segment_length);

    QVector<double> pressures(buffer.size()), temperatures(buffer.size());
    std::transform(buffer.begin(), buffer.end(), pressures.begin(), 
        [] (const VectorXd& v) -> double { return v(0); });
    std::transform(buffer.begin(), buffer.end(), temperatures.begin(), 
        [] (const VectorXd& v) -> double { return v(1); });

    plot->register_curve(0, 2e6, "pressure", Qt::red, QwtPlot::yLeft);
    plot->register_curve(50, 350, "temperature", Qt::blue, QwtPlot::yRight);

    plot->set_curve_data(QString("pressure"), &pipe.profile.coordinates.front(), &pressures.front(), buffer.size());
    plot->set_curve_data(QString("temperature"), &pipe.profile.coordinates.front(), &temperatures.front(), buffer.size());

    ui.Plot->replot();
	//graphicsDraw(ui.slTime->value());
}


void graph::on_slTime_valueChanged(){
	if (pde_time.size()){
		graphicsDraw(ui.slTime->value());
	}
	
}
void graph::graphicsDraw(int current_time){
    /*if (current_time >= pde_time.size()) {
        return;
    }

    QString s;
    s.setNum(pde_time[current_time]);
    ui.lbTime->setText(s);
    
    for (unsigned index = 0; index < 2; ++index) {
        if (current_time >= pde_data[index].size()) {
            continue;
        }
        list<vector<VectorXd> >::iterator U = pde_data[index].begin();
        std::advance(U, current_time);
        QVector<double> U0, U1, time;
        for (int i = 0; i < U->size(); i++){
            U0.append(U->at(i)(0));
            U1.append(U->at(i)(1));
            time.append(i);
        }
        curves[index][0]->setSamples(time, U0);
        curves[index][1]->setSamples(time, U1);
        if (index == 0)
            curves[0][2]->setSamples(&time[0],  &pipe_profile.heights.front(), time.size());

        if (index == 0) {
            ui.Plot->replot();
        }
    }
    for (unsigned index = 0; index < 2; ++index) {
        if (current_time >= scalar_data[index].size()) {
            continue;
        }
        list<vector<double> >::iterator U = scalar_data[index].begin();
        QVector<double> coordinates;
        for (int i = 0; i < U->size(); i++) 
            coordinates.append(i);

        std::advance(U, current_time);
        curves[0][index]->setSamples(&coordinates.first(), &U->front(), U->size());

        if (index == 0)
            curves[0][2]->setSamples(&coordinates.first(),  &pipe_profile.heights.front(), pipe_profile.heights.size());

        ui.Plot->replot();
    }*/
    
    
}
